Method for automatic branch labeling

ABSTRACT

The present invention relates to a method of analyzing an object data set in which a tubular structure having a plurality of branches and bifurcations occurs, wherein the object data set assigns data values to positions in a multi-dimensional space, which data values relate to an object to be examined. The invention relates further to a corresponding apparatus and computer program.

The present invention relates to a method of analysing an object dataset in which a tubular structure having a plurality of branches andbifurcations occurs, wherein said object data set assigns data values topositions in a multi-dimensional space, which data values relate to anobject to be examined. The invention relates further to a correspondingapparatus and computer program.

Nowadays, it is possible to acquire volume representations of the brainwhich show a clear distinction in grey values between tissue and vesselvoxels. Shape extraction, such as measuring a vessel's cross-sectionalarea, is done by interactively positioning and orienting a plane. Theintersection of this plane with the volume gives a 2D image of greyvalues in which the vessel pixels have a different grey value than thetissue pixels. After selection of the proper object the cross-sectionalarea can be measured, for example by counting the number of non-blackpixels. The plane should be oriented so that it is orthogonal to thevessel whose shape has to be measured. An oblique plane would give thewrong cross-sectional area. Unfortunately, interactively orienting theplane orthogonally to the vessel is a time-consuming and error-pronetask.

In J. Bruijns “Semi-Automatic Shape Extraction from Tube-like Geometry”,Proceedings Vision Modeling and Visualization 2000, Saarbruecken,Germany, pp. 347-355, November 2000 a new tool, the self-adjustingprobe, was introduced. A probe is a combination of a sphere and a planethrough the centre of the sphere. After the probe has interactively beenplaced on a vessel in the neighbourhood of the desired position, theprobe automatically adjusts itself so that its plane is orthogonal tothe vessel and the centre of its sphere is on the central axis of thevessel. The ellipse radii of the vessel are also estimated. When theprobe is aligned, it can be moved along the vessel in the direction ofthe plane normal. The probe aligns itself again after each step. It istherefore possible to let the probe follow the vessel automaticallyuntil the probe detects for example the end of a vessel or the beginningof an aneurysm.

Currently, the self-adjusting probe extracts shapes from a 3D trianglesurface model, created for example by a known marching cube algorithm.However, if two vessel branches are close together, it is possible thatsurface vertices of the neighbour vessel branch are included in the setof selected surface vertices which are used to extract the local shapeof the vessel branch investigated. These erroneously included surfacevertices deteriorate the accuracy of the estimated shape.

It is therefore the object of the present invention to provide a methodof analysing an object data set having an improved accuracy.

This object is achieved by a method comprising the steps of:

-   -   finding the extremities of the branches of said tubular        structure,    -   forming a skeleton of branches and bifurcations by a peeling        step,    -   forming directional graphs for the branches of said skeleton        between two neighbouring bifurcations or between a bifurcation        and an extremity based on said skeleton,    -   assigning a label to the positions along the directional graphs,        wherein for each branch of each directional graph a unique label        is selected,    -   determining the geometry of the branches and bifurcations of        said tubular structure so that positions can be classified as        belonging to either a bifurcation or a branch, and    -   assigning a final label to the positions along the branches and        of the bifurcations of said tubular structure, wherein for each        branch and each bifurcation a unique label is selected.

The invention is based on the idea to give the vessel voxels, and fromthese, the surface vertices, a unique label per vessel branch of saidtubular structure. Now, vertices of neighbour vessel branches can beexcluded because their label is different. This new method results alsoin a set of directed graphs—one for each component—which facilitatesfully automatic vessel tracing from one node, i.e. an extremity orbifurcation of the vessel structure, to another node of the same graph.

The fully-automatic branch labelling method according to the inventionconsists of five steps. Starting point is a segmented voxel volumewithout tissue inclusions. The elements of this segmented voxel volumeare for-example signed bytes, with a 0 for a tissue and a 1 for a vesselvoxel. This makes it possible to assign different labels to vesselvoxels during these steps. Vessel voxels with a label 1 are calledoriginal vessel voxels. Tissue voxels are never changed. The finaloutcome is a segmented voxel volume in which almost all vessel voxelshave a label indicating to which bifurcation or branch they belong and aset of directed graphs describing the topology of the vessels.

In a first step, the extremities of the voxel vessel structures aresearched, preferably by applying a wave propagation method, as e.g.described in Zahlten et al., “Reconstruction of Branching Blood VesselsFrom CT-Data”, Proceedings of the Eurographics Workshop on Visualizationin Scientific Computing, Rostock, June 1994. In a second step, thesegmented voxel volume with the extremities indicated by a special labelis peeled in a number of iterations. The resulting skeleton of branchesand bifurcations is a good approximation of the centre structure of thevessels. In a third step, the graphs are created. Starting point is thepeeled segmented voxel volume with the bifurcation and extremity voxelsindicated by special labels. A directed graph is generated for eachconnected set of vessel voxels in the peeled segmented voxel volume.This graph contains one node for each extremity and one node for eachbifurcation voxel. The vessel voxel strings between adjacent nodes arestored in branch lists. Each branch list gets a unique number. Thevoxels of a branch list get this number as label. In a fourth step, nodegeometry is generated. Node geometry contains the position and the shapeof the centre region of the bifurcation of said tubular structure andthe position and the shape of its branch regions. In a fifth step thevessel voxels get their final label. First, the voxels in the branchregions of said tubular structure get the label of the correspondingbranch list. Next, the voxels in the centre regions of said tubularstructure get a label which is different from the labels of the othercentre regions and different from the labels of the branch lists.Finally, the voxels of the branches between two adjacent branch regionsof said tubular structure get the label of the corresponding branchlist. After the vessel voxels are labelled, a modified marching cubealgorithm is applied which generates the surface vertices including thelabel of their neighbouring vessel voxel.

Preferred embodiments of the invention are defined in the specification.An apparatus for analysing an object data set is also defined in thespecification. The invention is further defined by a computer programfor implementing the method according to the invention.

According to a preferred embodiment the label of the correspondingbranch is assigned as final label to the positions along the branches.Mainly only for bifurcations new final labels have to be selected.

According to another embodiment the geometry of the branches and thebifurcations is determined on the basis of branch spheres andbifurcation spheres resulting in the position and shape of thebifurcations and adjacent branch regions. Thus, voxels in all regions ofbranches and bifurcations can correctly be labelled thus improvingaccuracy of the whole method.

Preferably, the branch spheres and bifurcation spheres for branches andadjacent bifurcations are arranged such that the radius of these spherescorresponds to the diameter of the corresponding branch or bifurcation,respectively. Further, the radius of a branch sphere or a bifurcationsphere is derived from the distance transform of the centre position ofthis sphere. Thus, voxels or bifurcation can be distinguished frombranch voxels which is required for correct labelling of voxels in asubsequent step.

The invention is generally applicable for analysis of amulti-dimensional object data set of a tubular object having a pluralityof branches and bifurcations, such as lung structures or vessel trees.The object data set can be acquired, for example, by means of varioustechniques such as magnetic resonance angiography, computed tomographyor 3D rotational X-ray angiography. Techniques of this kind produce anobject data set with data values representing the structure of (a partof) the vascular system of the patient to be examined.

The invention will now be explained in more detail with reference to thedrawings, in which

FIG. 1 shows a flow chart of the method according to the invention,

FIG. 2 shows a pixel map of a segmented voxel volume,

FIG. 3 illustrates wave propagation,

FIG. 4 illustrates splitting of waves,

FIG. 5 shows trial waves according to the invention,

FIG. 6 shows the segmented voxel volume with extremities indicated,

FIG. 7 shows a peeled segmented volume with extremity and bifurcationvoxels indicated,

FIG. 8 shows directional graphs obtained from the peeled segmentedvolume,

FIG. 9 illustrates the determination of note geometry,

FIG. 10 shows centre and branch spheres along the directional graphs,

FIG. 11 illustrates the labelling of voxels and

FIG. 12 shows a labelled segmented volume.

FIG. 1 shows a flow chart illustrating the steps of the method accordingto the invention. These steps will further be explained with referenceto the subsequent figures. In the first step S1 the wave propagationmethod is applied to a segmented voxel volume of a tubular structure, inthis particular example of a vessel structure of a patient, as exemplaryshown in FIG. 2. Conventionally wave propagation is used to label thegrey value voxel volume and to generate the corresponding vessel graph.However, wave propagation is not accurate enough at the bifurcationssince the incoming waves continue too long at the bifurcations.According to the invention wave propagation is used to generate theextremities of the vessel voxel structures. If the current wave is atthe far end of a branch, one of the voxels of the current wave isselected for inclusion in the set of extremities. The labels assignedduring the wave propagation step are ignored. The starting points aregenerated automatically. Selection by the user is too time-consumingbecause the segmented voxel volumes generally contain a lot of vesselgraphs. These starting points are included in the set of extremities.

A wave, with a branch number greater than 1, is a list of voxeldescriptions. A voxel description contains a label with the value thevoxel had when the description was created (except for trial wavesdiscussed later), the memory address of the voxel, giving fast access tothe current value of the voxel and its index (ix, iy, iz).

Starting with a single seed voxel as initial wave, a new wave with thesame branch number is created by generating voxel descriptions for allcorner neighbour vessel voxels, i.e. a voxel which shares at least onecorner with the voxel examined, of the current wave which were not yetmember of any wave. Whether a voxel was already member of a wave can bededuced from the current value of the voxel in the segmented voxelvolume because voxels which are included in a wave get the branch numberof the current wave. Next, the current wave is deleted and the new wavebecomes the current wave.

If the voxels of the current wave are not corner connected, the wave issplit in two or more new waves with each a new unique branch number.These new branch numbers are also assigned to the voxels of the splitwaves. Splitting will occur when a wave travels through a bifurcation ofthe vessel voxel structures. The wave propagation method is applied toeach of the new wave fragments. The process is finished when all newwaves are empty. The wave propagation method in general is explained indetail in the above mentioned article of Zahlten et al. In this articlethe generation of three waves is shown in FIG. 3 which is also includedherein as FIG. 3. FIG. 3A shows the initial wave consisting of a singleseed voxel (the dark cross). FIG. 3B shows the new wave consisting ofthe corner neighbour voxels (the five dark crosses) with the voxel ofthe initial wave marked by a filled circle. FIGS. 3C and 3D show thewaves generated from the corner neighbour voxels of the wave of FIG. 3Brespectively. FIG. 3D shows the splitting of a wave.

As already mentioned, wave propagation is preferably used—other methodare also possible—to find the extremities of the voxel structures.Extremities are found when the new wave is empty. However, in case of anoisy surface the wave will be split in many sub-waves just before theend of a branch is reached. This yields many unwanted extremities. Inthe example of FIG. 4 there are three waves shown. The first wavecontains the voxels labelled 1, the second wave the voxels labelled 2and the third wave the voxels labelled 3. The voxels of the third waveare not corner connected. So, this wave is split in three new waves.Because these three wave fragments do not have successors, threeextremities are generated.

Double waves solve this problem. The waves described in the previousparagraphs are called single waves in the following. A double wavecontains old voxels, already present in the previous double wave, andnew voxels, which are the corner neighbour vessel voxels of these oldvoxels. When from such a double wave a new double wave is created, thenew voxels of the old double wave are copied to the new double wave asold voxels. Next, the corner neighbour vessel voxels of the old voxelsof the new double wave are added as new voxels. The old voxels of a newdouble wave are also used to test whether the new voxels of this newdouble wave are corner connected. The new double wave is consideredempty if it contains only old voxels.

In the example of FIG. 4 there are three double waves. The first wavecontains the voxels labelled 1 as old and the voxels labelled 2 as newvoxels. The second wave contains the voxels labelled 2 as old and thevoxels labelled 3 as new voxels. The third wave contains the voxelslabelled 3 as old and no new voxels. The new voxels of the second doublewave are corner connected to each other via the old voxels of this wave.So, this second double wave is not split as was the case for the thirdsingle wave. The third double wave contains only old voxels. Hence, thiswave is considered empty and the extremity is derived from all voxels ofthe second double wave, giving only one extremity for this branch end.In the following sections, the waves for wave propagation (except thespecial waves) are double waves.

Wave propagation is based on corner neighbours. But peeling, graphgeneration, generation of node geometry and final labelling are based onface neighbours, i.e. a voxel which shares a face with the voxelexamined. Indeed, each voxel has 26 corner neighbours and 6 faceneighbours. So, using face neighbours instead of corner neighbours savesa significant amount of computing time. Even more important, cornerneighbours may yield much more unwanted connections between parallelneighbour vessels. Therefore, before wave propagation is started for anew seed voxel, the original vessel voxels which are face connected, viavessel voxels, to this seed voxel, are labelled as component voxels. Theseed voxel itself is also labelled as component voxel.

These face connected original vessel voxels are labelled by thefollowing twin wave algorithm: Two empty waves are created. The firstwave is filled with a voxel description of the seed voxel. The seedvoxel is labelled as a component voxel. The second wave is filled withvoxel descriptions of those original vessel voxels which are faceneighbours of the voxels of the first wave. The corresponding voxels arelabelled as component voxels. This process is repeated, changing eachtime the role of the first and second wave, until no original vesselvoxels, face connected to the voxels of the current wave are found. Thenormal wave propagation can now be applied on the component voxels.Vessel voxels which are corner but not face connected to a voxel of thiscomponent are now ignored because they do not have the component label.

Wave propagation yields the best results, if the waves travel from widerto narrower vessels. Because wave propagation is used to find theextremities of the vessel structures, seed voxels should be in theneighbourhood of these extremities. Fortunately, very often the largedominant vessel structures start on the boundaries of the volume withtheir widest vessels. The centre voxels of the widest vessel can befound by means of the distance transform. The distance transform of avessel voxel is an indicator for the length of the shortest path of faceconnected vessel voxels to a voxel of a vessel boundary.

It is possible that a local maximum in the distance transform is sharedby a number of face neighbour vessel voxels in the same boundary volumeslice. So, it is not safe to look for the vessel voxel with a distancetransform greater than the distance transform of its face neighbours.Therefore, all voxels in the six boundary volume slices with a distancetransform greater than or equal to the distance transform of its faceneighbours are stored in a special wave, the so-called seed wave, bymeans of an insertion sort on the basis of their distance transforms.The insertion sort makes that the seed voxel with the highest distancetransform is the first voxel description in the seed wave.

When a voxel is included in the seed wave, the distance transform of itsface neighbour vessel voxels is decreased. This prevents that manyvessel voxels with the same distance transform as their face neighboursare included in the seed wave. Starting with (and removing) the firstvoxel description of the seed wave, the wave propagation algorithmincluding the component selection is carried out. Many voxels of theseed wave will be altered because they belong to the same vesselcomponent. So, skipping and removing the first voxel descriptions of theseed wave from which the corresponding voxels have already been changed,the wave propagation algorithm including the component selection iscarried out for all first voxel descriptions of the seed wave from whichthe corresponding voxels are not yet processed, until the seed wave isempty.

Many components do not start or end on a boundary of the volume.Scanning the segmented voxel volume for not yet processed voxels, willnot always result in a seed voxel in the neighbourhood of theextremities of a component. But, by sending trial waves from the startvoxel found by scanning, the extremities will be found as shown in FIG.5. Therein the volume boundary is indicated by 5, the volume scandirection by 6, the left trial path by 7 and the right trial path by 8.The extremities 4 of the trial wave propagation pass are stored in theseed wave instead of in the set of extremities.

Waves do not only have a branch number but also a serial number. Thisserial number is increased every time a new wave is created from an oldwave. The voxel descriptions of the extremity voxels of the trial wavesget the serial number of the trial wave as label (not the current valueof the extremity voxel as is normally the case). Because the seed waveis filled by means of an insertion sort, the extremity voxel of thetrial wave with the highest serial number will be used as seed voxel fora normal wave propagation pass. The trial waves get a special branchnumber, not used in case of a normal wave propagation pass. So, afterthe seed voxels are found, the voxels of the component can be reset totheir original value. The trial wave propagation pass is alternated withthe normal wave propagation pass until all interior components have beenprocessed.

The result of the wave propagation step S1 is a segmented volume withextremity voxels as shown in FIG. 6 where the extremities are indicatedas darker spots 9.

In the second step S2, the segmented voxel volume is peeled in a numberof iterations by an algorithm similar to that described by Dokládal etal. “A new thinning algorithm and its application to extraction of bloodvessels”, conference proceedings, BioMedSim '99, ESIEE, April 1999,France. The resulting skeleton of branches and bifurcations shown inFIG. 7 is a faithful approximation of the centre structure of thevessels: the branches are close to the core lines of the vessel branchesand the skeleton bifurcations are close to the centres of the vesselbifurcations.

There are many algorithms to create such a skeleton. In the followingone particular algorithm shall be explained, which may, however, bereplaced by another algorithm as long as the resulting skeleton is afaithful approximation of the centre structure of the vessels.

Starting point is the original segmented voxel volume with the extremityvoxels indicated by a special label, the darker spots 9 in FIG. 6. Eachiteration creates first a skin layer by labelling the current boundaryvoxels except the extremity voxels. Each voxel has at most six faceneighbours. An original vessel voxel with label one is a boundary voxelif and only if one of these neighbours is not a voxel with a positivelabel. The label of a boundary is increased with the number of faceneighbours with a zero or negative label (boundary voxels removed in aprevious iteration).

Checking and removing (if the local topology does not change) allboundary voxels of the current skin layer before creating, checking andremoving the boundary voxels of the next inner skin layer, guaranteesthat the remaining set of voxels approximates the core lines of thevessel graphs. After peeling is finished, bifurcation voxels are markedby a special label. Bifurcation voxels are voxels with more than twopositive neighbours in the peeled segmented voxel volume shown in FIG.7.

The boundary voxels of the current skin layer are checked and possiblyremoved in order of their label. The volume is traversed from the lowerleft front corner to the upper right back corner, looking for boundaryvoxels with the current label. Boundary voxels with exactly one vesselvoxel as face neighbour are always removed. Boundary voxels with atleast two and at most four vessel voxels as face neighbours are removedunless the local topology changes by removing this voxel. Boundaryvoxels with exactly five vessel voxels as face neighbour are neverremoved because removing such a boundary voxel causes a localcon-cavity. If a boundary voxel is removed, the labels of its faceneighbour boundary voxels are adjusted. Face neighbours with a labelgreater than or equal to the label of the boundary voxel removed, areimmediately processed in decreasing order of their new label beforenormal scanning from the lower left front corner to the upper right backcorner continues.

The topology checks are performed by looking at the corner neighbours ofthe boundary voxel tested. To this end a cube of 3 by 3 by 3 cells isfilled with zero's. The centre cell corresponds to the boundary voxelchecked. The other cells correspond to its corner neighbours. A cell isset to one if the corresponding corner neighbour has a positive label.It is possible that some of these cells are equal to one although thecorresponding corner neighbour vessel voxel is not face connected, viavessel voxels, to the boundary voxel examined. This can happen when twocomponents are corner but not face connected to each other. These cellsare reset to zero. First, resetting the centre cell to zero, should notresult in two disjunct face connected sets of positive cells. Secondly,resetting the centre cell to zero, should not decrease the number offace or the number of corner connected sets of zero cells. If one ofthese two complementary checks fails, the boundary voxel examined shouldnot be removed.

In the third step S3, the graphs for the branches of said skeleton arecreated. Starting point is the peeled segmented voxel volume as shown inFIG. 7 with the bifurcation voxels indicated by a special label togetherwith the extremities, stored in a special wave by the wave propagationstep. A directed graph is generated for each face connected set ofpositive vessel voxels in the peeled segmented voxel volume. This graphcontains one node for each extremity voxel and one node for eachbifurcation voxel. A wave is created and stored in this graph for eachbranch between two nodes. This special wave consists of a list, i.e. abranch list, of voxel descriptions for the face connected positivevessel voxels between these two nodes. Each branch list gets a uniquenumber. The voxels of a branch wave (branch list) get this number asspecial unique label. An example of the resulting segmented voxel volumeis shown in FIG. 8.

The generated graphs facilitate not only fully-automatic vessel tracingfrom one node to another node, but are also required to label shortbranches correctly. In this last case, information about the bifurcationstructure, especially its size, at one end of the branch is needed forlabelling voxels at the other end of the branch. The generated graphscontain this neighbour information.

A graph contains amongst other things the number of branches in thegraph, the number of nodes in the graph, a pointer to the list of itsbranch waves and two pointers, one to the first node and one to the lastnode of the list of nodes of the graph. A node contains amongst otherthings a voxel description for the corresponding extremity orbifurcation voxel (called the centre voxel of the node in the following)and the number of branches in this node. The number of branches shouldbe either equal to one or greater than two and less than or equal tosix. If the number of branches is equal to one, the node corresponds toan extremity of the vessel structure. Two branches at a node is notpossible because nodes are created either for an extremity or abifurcation. A bifurcation voxel is a voxel with at least three and atmost six positive face neighbours.

A node contains for each branch amongst other things a pointer to thecorresponding branch wave, a direction number indicating whether thehead or the tail of this branch wave is connected to this node and apointer to the node at the other end of this branch wave. The pointersto other nodes represent the graph structure. The branch waves make itpossible to travel from one node to another node going from one voxel tothe next voxel of the branch wave.

It should be noted that whether a branch of a node is incoming oroutgoing does not imply direction of blood flow at that node.

The graphs are generated by scanning the list of extremities. It is easyto detect whether an extremity is already part of a graph because duringthe generation of the graph, the voxels visited are made negative. Foreach positive extremity voxel a node is generated and filled. For eachpositive face neighbour of the centre voxel of a node a branch wave isgenerated and filled with the face connected positive vessel voxelsuntil a positive extremity or bifurcation voxel is encountered. A nodeis generated for the closing extremity or bifurcation voxel and thebranch information in the two nodes is updated. The algorithm continueswith tracing the branches of the new node.

In case of cycles in the vessel structures, a node may already exist forthe closing bifurcation voxel. If so, the bifurcation voxel is negative.Therefore, if no positive extremity or bifurcation voxel is found at theend of a branch wave, the face neighbours of the last positive voxel areinspected for a negative bifurcation voxel. The corresponding node isused to close the branch.

Correct labelling of voxels requires that voxels of a bifurcation can bedistinguished from branch voxels. To this end node geometry of thebranches and bifurcations of the tubular structure is generated in stepS4 as shown in FIG. 9 so that positions can be classified as belongingto either a bifurcation or a branch. In this step the vessel boundary isindicated by 18. Node geometry contains the position and the shape ofthe centre region of the bifurcation and the position and the shape ofits branch regions. First, a centre sphere 10 is created. The positionof the bifurcation voxel is used as centre 11 of the centre sphere 10.The radius (in voxels) of the centre sphere is derived from the distancetransform of the bifurcation voxel.

Next, a branch sphere 12 is created for each branch 13. The centre 14 ofa branch sphere 12 is equal to the position of the voxel of the branchwave so that the branch sphere 12 is just separated from the centresphere 10. The radius of the branch sphere 12 is derived from thedistance transform of the centre voxel of the branch sphere 12.Travelling along and checking each voxel of the branch wave yields thefirst voxel which fulfils these conditions.

Finally, one branch plane 15 and one centre plane 16 is created for eachbranch 13. The branch plane 15 is defined by the centre 14 of the branchsphere 12 and the normal 17 which is given by the direction of theconnection line from the centre 14 of the branch sphere 12 to the centre11 of the centre sphere 10. The position of the corresponding centreplane 16 is determined by the intersection of the centre sphere 10 andthis connection line 17. Its normal is equal to the normal of the branchplane 15 multiplied with −1.

If the distance in voxels along the connecting branch wave between twobifurcation voxels or between a bifurcation and extremity voxel is smallcompared with the radii of the centre spheres, centre sphere 10 andbranch spheres 12 overlap. Identifying the entities at one end with“first” and at the other end with “second”, the following cases can bedistinguished:

-   1. If all voxels of a branch wave are inside the first centre sphere    10, the position of the last voxel is used as centre of the first    branch sphere 12. The radius of this branch sphere 12 is in this    case multiplied with minus one to indicate this condition.-   2. If some of the voxels of this branch wave are outside the first    centre sphere 10, but the first branch sphere 12 overlaps the first    centre sphere 10 even for the last voxel of the branch wave, this    last voxel is used as centre of the first branch sphere 12. This    overlap does not harm the final labelling because only the voxels    between the branch plane 15 and the corresponding centre plane 16    are members of the branch region.-   3. If the second node of the branch is a bifurcation and if the    centre of the first branch sphere 12 is inside the second centre    sphere, the radius of the first branch sphere 12 is multiplied with    minus one to indicate this condition. Indeed, if the centre of a    branch sphere 12 is inside one of the centre spheres 10 of a branch    between two bifurcations, all voxels of the branch 13 are inside the    two bifurcations.

It should be noted that if the node at the other end of the branch 13 isan extremity, all voxels outside the bifurcation including the voxels inthe neighbourhood of the extremity are labelled as branch voxels. So, inthis case it is not necessary to test whether the centre of the branchsphere 12 of the bifurcation is inside the centre sphere of theextremity.

The branch and centre spheres are shown in FIG. 10. As is clear fromthis figure, centre spheres are also generated for the extremities. Theoccasional overlap of spheres is also shown in this figure.

In the fifth step S5, the vessel voxels get their final label. Startingpoint is the original segmented voxel volume with only tissue andoriginal vessel voxels and the node geometry generated in step S4. Theoriginal vessel voxels in a region are labelled by a similar twin wavealgorithm as is used for component selection as described above. Thedifference is that the original vessel voxels should not only be faceconnected, via vessel voxels, to the initial voxel, but should alsofulfil constraints depending on the region currently processed.

First, the voxels in the branch regions 19 of the tubular structure getthe label of their branch 13, i.e. of the corresponding branch list. Asimplified 2D example of a branch region 13 is given in FIG. 11. Theinitial voxel for the twin wave algorithm is the centre voxel 14 of thebranch sphere 12. The additional constraints are:

-   1. The distance of the original vessel voxel to the centre 14 of the    branch sphere 12 should be less than or equal to twice the radius of    the branch sphere 12.-   2. The original vessel voxel should reside between the branch plane    15 and the centre plane 16. The first constraint prevents    unrestricted growth of the labelled area of a side branch 13 in case    the centre sphere 10 is so small that the centre plane 16 for this    side branch 13 does not intersect this branch 13 (as to the left of    the right vertical vessel boundaries in FIG. 11). In this case a    large number of original vessel voxels of the main branch 13′,    fulfilling the second constraint, are face connected to the centre    voxel 14 of the branch sphere 12.

After the voxels in the branch regions 19 are labelled, the voxels inthe centre regions 10 of the tubular structure are labelled. Each centreregion 10 gets a unique label, different from all branch numbers anddifferent from the labels assigned to the other centre regions and thebranch lists. The initial voxel for the twin wave algorithm is thecentre voxel 11 of the node. The additional constraints are:

-   1. The distance of the original vessel voxel to the centre 11 of the    centre sphere 10 should be less than or equal to the radius of the    centre sphere 10 plus the maximum of the branch radii of the current    node.-   2. The original vessel voxel should reside inside the “enclosure” of    the branch planes 15 of the node.-   3. The distance of the original vessel voxel to the position of the    current node should be less than or equal to all distances of the    original vessel voxel to the positions of the neighbour nodes.

After labelling the original vessel voxels in the branch regions 19, theoriginal vessel voxels of the centre regions 10 are separated from theremaining original vessel voxels in the branches 13 by the alreadylabelled vessel voxels. But when a branch region 19 is skipped becauseits branch radius is negative or when the distance between the branchplane 15 and centre plane 16 is very small, the original vessel voxelsof the centre region 10 are face connected to the original vessel voxelsof the branch 13. The first two constraints prevent unrestricted growthof the labelled area in these cases. The third constraint separates thevoxels of the current centre region from the voxels of the centreregions of neighbour nodes in case of overlapping centre regions.

After the voxels in the centre regions 10 are labelled, the voxels ofthe branches 13 between two branch regions 19 of the tubular structureare labelled and get the label of the corresponding branch list. Theoriginal vessel voxels of the branch wave are used as initial voxels.But using one original vessel voxel of the branch wave as initial voxelfor the twin algorithm will normally result in labelling of most otheroriginal vessel voxels of the branch wave. So, only one or two vesselvoxels of the branch wave will really be used as initial voxel.

Because the remaining original vessel voxels of a branch are separatedfrom the original vessel voxels in other branches by the labelled voxelsof the centre regions, no additional constraints are necessary toprevent unrestricted growth of the labelled area. A labelled segmentedvolume is shown in FIG. 12.

The following conclusions can be drawn when applying the methodaccording to the invention:

-   1. The new method for fully-automatic branch labelling of voxels in    vessel structures gives better results than the wave propagation    method.-   2. The wave propagation step deteriorates in case of wide short    branches (for example an aneurysm), producing spurious extremities.-   3. The quality depends on the smoothness of the surfaces of the    voxel structures. A very noisy surface results in many short    branches. Smoothing the voxel volume can remove these short    branches. However, finding the correct smoothing factor may still    require human interaction.-   4. The elapsed time for the wave propagation steps depends on the    number of components. The elapsed time for the peeling step depends    on the number of vessel voxels in the original segmented volume. The    elapsed time for the last three steps is negligible compared with    the time taken by the first two steps.-   5. The method according to the invention results also in a set of    directed graphs, one for each component, which facilitates    fully-automatic vessel tracing from one node—an extremity or    bifurcation of the vessel structure—to another node of the same    graph.

1. A method of analyzing an object data set in which a tubular structurehaving a plurality of branches and bifurcations occurs, wherein saidobject data set assigns data values to positions in a multi-dimensionalspace, which data values relate to an object to be examined, said methodcomprising the steps of: finding the extremities of the branches of saidtubular structure, forming a skeleton of branches and bifurcations by apeeling step, forming directional graphs for the branches of saidskeleton between two neighboring bifurcations or between a bifurcationand an extremity based on said skeleton, assigning a label to thepositions along the directional graphs, wherein for each branch of eachdirectional graph a unique label is selected, determining the geometryof the branches and bifurcations of said tubular structure so thatpositions can be classified as belonging to either a bifurcation or abranch, and assigning a final label to the positions along the branchesand of the bifurcations of said tubular structure, wherein for eachbranch and each bifurcation a unique label is selected, and wherein thegeometry of the branches and the bifurcations is determined on the basisof branch spheres and bifurcation spheres resulting in the position andshape of the bifurcations and adjacent branch regions.
 2. The method ofclaim 1, wherein the branch spheres and bifurcation spheres for branchesand adjacent bifurcations are arranged such that the radius of thespheres corresponds to the diameter of the corresponding branch orbifurcation, respectively.
 3. The method of claim 1, wherein the radiusof a branch sphere or a bifurcation sphere is derived from the distancetransform of the centre position of the sphere.